#include "solver3.h"
#include <math.h>

Solver3::Solver3(double _u0, double _du0, double _a, double _b, double r, double eps, double eps_r, int n, double _h): u0(_u0), du0(_du0), b(_b), epsilon(eps), epsilon_right(eps_r), Nmax(n), h(_h), right(r), a(_a)
{
	result1 = new std::list<SearchNode>();
	result2 = new std::list<SearchNode>();
}

double Solver3::f1(double x, double v, double u)
{
	return u;
}

double Solver3::f2(double x, double v, double u)
{
	return (-1)*(a*a*sin(v)+b*sin(x));
}

void Solver3::solve()
{
	result1->push_back(SearchNode(h,  0.0, u0, 0.0, 0.0, 0.0, u0, u0));
	result2->push_back(SearchNode(h,  0.0, du0, 0.0, 0.0, 0.0, du0, du0));

	double x = 0.0;
	double v = u0;
	double u = du0;
	
	int n = 0;

	while ((x < (right - epsilon_right)) && (n < Nmax))
	{
		double x_new = x + h;
		double v_new = v + h/2.0*(f1(x,v,u) + f1(x+h,v+h*f1(x,v,u), u+h*f1(x,v,u)));
		double u_new = u + h/2.0*(f2(x,v,u) + f2(x+h,v+h*f2(x,v,u), u+h*f2(x,v,u)));

		double v_mid = v + h/4.0*(f1(x,v,u) + f1(x+h/2.0,v+h/2.0*f1(x,v,u), u+h/2.0*f1(x,v,u)));
		double u_mid = u + h/4.0*(f2(x,v,u) + f2(x+h/2.0,v+h/2.0*f2(x,v,u), u+h/2.0*f2(x,v,u)));
		
		double v2 = v_mid + h/4.0*(f1(x + h/2,v_mid, u_mid) + f1(x+h,v_mid+h/2*f1(x + h/2,v_mid, u_mid), u_mid+h/2*f1(x + h/2, v_mid, u_mid)));		
		double u2 = u_mid + h/4.0*(f2(x + h/2,v_mid, u_mid) + f2(x+h,v_mid+h/2*f1(x + h/2,v_mid, u_mid), u_mid+h/2*f1(x + h/2, v_mid, u_mid)));
		
		double s1 = (v2 - v_new)*4/3;
		double s2 = (u2 - u_new)*4/3;

		double s = std::max(abs(s1),abs(s2));

		if (abs(s) < epsilon) {
			if (abs(s) > epsilon/8.0) {
				result1->push_back(SearchNode(h, x_new, v_new, v_mid, v2, s1, v_new, 0));
				result2->push_back(SearchNode(h, x_new, u_new, u_mid, u2, s2, u_new, 0));
			}
			else {
				result1->push_back(SearchNode(h, x_new, v_new, v_mid, v2, s1, v_new, 0));
				result2->push_back(SearchNode(h, x_new, u_new, u_mid, u2, s2, u_new, 0));
				h = h*2;
			}
			x = x_new;
			v = v_new;
			u = u_new;
		}
		else {
			h = h/2;
		}
		n++;
	}
}
void Solver3::getPortret()
{
	result1->clear();
	result2->clear();
	
	double x = 0.0;
	double v = u0;
	double u = du0;

	for (x = 0.0; x < right; x+=h)
	{
		double v_new = v + h/2.0*(f1(x,v,u) + f1(x+h,v+h*f1(x,v,u), u+h*f1(x,v,u)));
		double u_new = u + h/2.0*(f2(x,v,u) + f2(x+h,v+h*f2(x,v,u), u+h*f2(x,v,u)));

		result1->push_back(SearchNode(0.0, v_new, 0, 0, 0, 0, u_new, 0));
		v = v_new;
		u = u_new;
	}

	v = u0;
	u = du0;
	//for (x = 0.0; x > -right; x-=h)
	//{
		//double v_new = v - h/2.0*(f1(x,v,u) + f1(x-h,v-h*f1(x,v,u), u-h*f1(x,v,u)));
		//double u_new = u - h/2.0*(f2(x,v,u) + f2(x-h,v-h*f2(x,v,u), u-h*f2(x,v,u)));

//		result2->push_back(SearchNode(0.0, v_new, 0, 0, 0, 0, u_new, 0));
	//	v = v_new;
		//u = u_new;
	//}
}